N dimensional non-linear, static, adaptive digital filter design using d scale non-uniform sampling

ABSTRACT

The present invention is a filter design that extracts information from a signal by employing D scale nonuniform sampling. In one embodiment a D scale multiresolution sampler, filter bank router, filter bank sampler controller, phase shifter, and consolidator constitute a maximal arrangement for a D scale FIR/IIR filter design.

This application is a Continuation of U.S. patent application Ser. No.10/250,830, filed on Jan. 30, 2004, which is now U.S. Pat. No. ______issued on ______, which is the National Phase of InternationalApplication No. PCT/US02/00940 filed Jan. 4, 2002, which is based onApplication Ser. No. 60/259,961, filed Jan. 5, 2001.

FIELD OF THE INVENTION

This invention concerns N dimensional digital filter design in itsbroadest sense. Filters are used to restrict or sculpt 1D signals, 2D/3Dimages or more abstract N dimensional objects that traversecommunications channels. They are also used for data smoothing andprediction, signal detection and noise reduction applications.

Typically, digital filters rely on the Finite Impulse Response (FIR) andInfinite Impulse Response (IIR) linear mathematical models. Data inputsto these models are almost always assumed to be uniformly spaced samplepoints. Various filter architectures use these underlying models toaccommodate diverse application categories, to include:

a. Static filtering (e.g., baseband, passband, notch)

b. Adaptive and equalizing filters

c. Weiner or matched filtering

d. Kalman filtering

This invention also concerns novel uses of digital filters to modulatepositive and negative feedback loops in control systems as well as toenhance digital phase lock loops.

Further, this invention extends the domain of applicability of digitalfilter techniques to a certain range of nonlinear systems.

Finally, this invention encompasses the derivative work of exactlycalculating the n-th level autocorrelation and cross-correlations of Ndimensional signals as well as convolutions.

BACKGROUND OF THE INVENTION

This discussion is confined to one dimensional (1D) filters because theN dimensional case is readily generalized from the 1D case.

A 1D filter is characterized by a vector of possibly complex valuedweights or parameters. These weights are the tap values in a transversalfilter arrangement. These weights are derived from the discreteapproximation to the continuous convolution of an input signal and theimpulse response of the system to filter. The convolution characterizesthe filter's impact on the input signal, as follows:

y(t) = ∫_(−∞)^(+∞)h(t − τ)u(τ)τ

where u(τ) is the input signal and h(t−τ) is the (non-causal in thiscase) impulse response of the system shown in FIG. 1.

The impulse response of a system is the result of applying a sharpspiked signal at the input. The spiked signal is a discreteapproximation to the unrealizable impulse of zero duration and infiniteamplitude (defined using the Dirac delta functional). The idealizedspectrum of impulse is flat with unit amplitude and constant phase.

The discrete approximation of the impulse response is:

${y(t)} = {{\int_{- \infty}^{+ \infty}{{h\left( {t - \tau} \right)}{u(\tau)}{\tau}}} \approx {{1/N}{\sum\limits_{k = 0}^{k = {N - 1}}{{h\left( {n - k} \right)}{u(k)}}}}}$

where 1/N is the uniform, Nyquist, sampling rate used to digitize thecontinuous signal.

Digital filters for linear systems are categorized as Finite-ImpulseResponse (FIR) or Infinite Impulse Response (IIR) in architecture. TheIIR occurs when the output at one sample time is fed back as input tothe next sample event. FIGS. 2 and 3 depict a general purpose FIR or IIRdigital filter, respectively, with both architected as a transversefilter.

There are many varied techniques for designing digital filters. Staticfilter design entails the determination of static tap weights forspecific problem domains, bounded by requirements and resourceconstraints. Adaptive digital filters result when the tap weights areallowed to vary, typically depending on a comparison of the actualoutput to a reference output. Typically, the Weiner-Hopf equationprovides the analytical basis for calculating an adaptive filter'soptimal tap weights.

CURRENT STATE-OF-THE-ART Description of Known Solutions

The current state of the art in digital filter design encompasses a widerange of solutions to the underlying equations that determine the filterweights. These techniques are generally categorized by the filterarchitecture, to include:

a. Linear Filter Design Algorithms

b. Somewhat Nonlinear Filter Design Algorithms

c. Adaptive and Equalizing Filters

d. Weiner or Matched Filters

e. Kalman Filters

DEFICIENCIES OF EXISTING “PRIOR ART”

The prior art is generally based on the following premises:

a. Linear systems models.

b. Approximation of the Fourier Transform via the Discrete FourierTransform (DFT). This is a fair approximation under best signalcondition, and potentially considerably worse for transient andnon-stationary signals.

c. Generally periodic, stationary signals.

d. Decoupling of sampling from filtering operations. That is, samplingis performed with whatever approach independently of the filteringresults.

e. Uniform sampling of the input signal.

BRIEF DESCRIPTION OF THE INVENTION

The N dimensional digital filter design is generalized in astraight-forward manner from the 1D case. Therefore the description willfocus on the 1D case. An N dimensional object will then be referred toas a signal, using the 1D terminology.

The motivation for this invention evolved from a realization that acritical component was excluded from the typical information flowthrough digital filters. This flow is illustrated in FIG. 4. Inputs cantypically include the current input sample, previous samples, a prioriknowledge and reference data. But the A/D sampler is not directlyincorporated into the actual filter. Thus, standard filter designtypically begins with a vector of uniformly sampled input data.

Fundamentally, non-uniform sampling of a volatile signal can glean moreinformation about a signal than uniform sampling at the Nyquist rate (ifindeed that rate is known a-priori). Additional information might begleaned from simultaneous pseudo-random non-uniform sampling, along withtrend or volatility sampling. Relative to judicious non-uniformsampling, uniform sampling can be sub-optimal in collecting informationper unit of effort expended.

This filter design invention, named the D digital filter, extracts thatadditional information from the non-uniform sampling of signalvolatility to improve the effectiveness and performance of the filter'sintended application. The performance measure or metric used depends onthe application. For example, for standard noise filtering and signalsmoothing applications, the D filter is designed to sharpen spectrumcutoff and to minimize ripples in the frequency spectrum, whilepreserving phase. Or, for adaptive applications, the D filter isdesigned to accurately or closely track changes from a referencecondition, with fast response time. FIG. 5 illustrates the informationflow through a filter with sampling incorporated.

This invention is of a digital filter design comprised of up to seven(7) components, shown in FIG. 6. The components actually used in anyimplementation will depend on the application. Typically, adaptivefiltering design will use many of the design components while staticfilters the fewest number.

The D filter design components are:

A. D Scale Multi-Resolution Input Sampler

B. Filter Bank Router

-   -   Filters can be applied to output streams following traversal        through a channel as well as to inputs. Each point in the output        stream might be shifted relative to its position in the input        stream. This component then maps each point to a transversal        filter in the bank of filters that is associated with a        particular D subscale.

C. Digital Filter Bank

-   -   Each filter in the bank implements a vector of tap weights using        the standard FIR and IIR architectures shown in FIGS. 2 and 3        respectively. Each filter is associated with one D subscale, and        the number of taps is the size of the D subscale. Collectively,        the filter bank then implements the D Scale underlying the        sampling.

D. Consolidator

-   -   This component consolidates or interleaves the multiple,        uniformly sampled D Subscale streams into one non-uniformly        spaced stream.

E. Arithmetic Mappings (Networked Tables)

-   -   This optional component is used to implement the various real        time calculations on the sample values.

F. Sampler Controller

-   -   This component sends control signals back to the sampler, to        dynamically change the sampling behavior on the input.

G. Phase Shifter

-   -   This component modifies the phase of output samples for each D        subscale (via a delay). A motivation is to modulate or control        positive, negative and phase lock feedback loops.

Several illustrations are presented to highlight the scope ofmix-and-match of these components, FIG. 7 shows the static D filterdesign. And FIG. 8 shows the Adaptive D filter design, both derived fromthe D filter components.

Further, the D Filter components might be used to modulate the inputsample streams based on the output. This can be implemented usingcontrol signals that are based on the output, and then applied to theinput samplers. The applied control signals would then suppress ortrigger certain D subscale samplers in various frequency ranges. TheSampler Controller module can serve that function. FIG. 8 shows a Dfilter arrangement that implements this sampling feedback loop.

Also, it might be desirable to modify the phase of an output or an inputsignal. The Phase Shifter component could fulfill that purpose, as shownin FIG. 9.

A D Filter can be applied in many configurations, identically to thestandard digital filters. FIG. 11 shows several possible static andadaptive D filter configurations.

The D filter design derives directly from the capabilities of the DScale (patent application Ser. No. 09/326,084), non-uniform samplingusing the D Scale, the D Arithmetic (patent application Ser. No.09/568,368) and on the exact Fourier transform decomposition of a signalvia D Scale sampling. Collectively, they provide the rigorousmathematical modeling justification for the D filter structure andoperation.

A signal structure is also introduced to aid the characterization of achannel for the construction of a D digital filter. Further, industrystandard static and adaptive digital filter designs are special cases ofthe D filter design when only uniform sampling is used. For then, thesampler is simply an A/D device, the filter bank reduces to one filterand there is no routing to a filter bank.

The D filter is also an evolutionary design, for it extends andextrapolates current filter design techniques. This preserves theinvestment in the vast installed base of digital filters.

Highlight of the D digital design include:

1. Enhanced IIR digital filtering by increasing its stability whileminimizing the number of tap weights. Stability is the main constraintto wider applications of current state-of-the-art IIR filters.

2. Sharper or finer “sculpting” of a channel's bandwidth to achieve thesame filtering effect as the usual designs, but with reduced number oftap weights. Alternately, the same number of taps as in a standardfilter yields heightened filtering impact.

3. Capability to retrofit existing filter designs because of itsbackward compatibility.

4. Evolutionary design is built on the existing framework of currentdigital filters. This reduces implementation risk as well as to preservethe investment in the installed base.

5. Mix-and-match component approach lends itself to solving many variedproblem in diverse application domains.

6. D Arithmetic self-stabilizing arithmetic further enhances the

7. stability of IIRs and of adaptive filters.

8. Efficiencies due to fine level of input sample filtering as well asfeedback from output back to input sampler.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 depicts an Industry Standard Linear System Model.

FIG. 2 depicts a Traversal FIR Digital Filter.

FIG. 3 depicts a Traversal IIR Digital Filter.

FIG. 4 depicts general information flow through a digital filter.

FIG. 5 depicts the general information flow through a filterincorporating a D scale sampler.

FIG. 6 depicts a maximal arrangement of D FIRR/IIR filter designcomponents.

FIG. 7 depicts a static D filter design.

FIG. 8 depicts an adaptive D filter design.

FIG. 9 depicts a D filter with a sampler controller.

FIG. 10 depicts a D filter with a phase modulator for positive, negativeand phase lock feedback loops.

FIG. 11 depicts possible D filter configurations.

FIG. 12 depicts multi-resolution D scale sampler designs.

FIG. 13 depicts a D digital filter band and sample router.

FIG. 14 depicts a filter bank router.

FIG. 15 depicts an internal architecture of an FIR D filter bank.

FIG. 16 depicts a D signal description.

FIG. 17 depicts the characterizing of certain nonlinear systems for Ddigital filter design.

FIG. 18 depicts a frequency domain analysis of nonlinear systems.

FIG. 19 depicts a procedure to characterize nonlinearities ininput-output amplitude or phases.

FIG. 20 depicts a methodology for empirically measuring the impulseresponse of a linear or nonlinear system.

FIG. 21 depicts a standard digital filtering model for Weiner-Hopfderivation.

FIG. 22 depicts a standard geometric interpretation of a solution to theWeiner-Hopf equation.

FIG. 23 depicts a digital filtering model for a generalized Weiner-Hopfderivation.

FIG. 23A depicts a standard continuous Weiner filter architecture.

FIG. 24 depicts a geometric interpretation of a solution to thegeneralized Weiner-Hopf equation.

FIG. 25 depicts a generalized Weiner matched filter.

FIG. 26 depicts a signal flow graph underlying a standard linear Kalmanfilter derivation.

FIG. 27 depicts a signal flow graph underlying a generalized Kalmanfilter derivation.

FIG. 28 depicts a generalized Kalman filter loop.

DETAILED DESCRIPTION Contents

1. Detailed Description of D Digital Filter Components

-   -   1.1 D Scale Multi-Resolution Input Sampler    -   1.2 Digital Filter Bank    -   1.3 Sample Router    -   1.4 Multi-Stream Consolidator    -   1.5 Phase Shifter    -   1.6 Sampling Feedback Loop    -   1.7 D Arithmetic Tables

2. Rationale for D Digital Filter

-   -   2.1 Decomposition of Convolution Integral    -   2.2 Internal Architecture of D Filter Bank    -   2.3 D Signal    -   2.4 Nonlinear Systems Addressed by D Filter Design    -   2.5 Method for Empirically Measuring Impulse Response of a        Linear or Nonlinear System

3. Application of D Filter Design

-   -   3.1 Linear Filter Design Algorithms    -   3.2 Nonlinear Filter Design Algorithms    -   3.3 Filters Derived from Weiner-Hopf Equation    -   3.4 Benefits of Generalized Weiner-Hopf Equation on Adaptive        Filters    -   3.5 Computing Autocorrelation and Cross-Correlation    -   3.6 Generalized Kalman Filters

1. Detailed Description of D Digital Filter Components

1.1 D Scale Multi-Resolution Input Sampler

Patent application Ser. No. 09/326,084 describes how the D Scale samplercan be implemented by multiple off-the-shelf, commercially availableanalog-to-digital (A/D) converters. Each A/D device samples uniformly ata rate determined by that D Scale's subscale resolution. An extracomponent is required to distribute or dispatch or fan out the originalanalog data stream into multiple feeds, one per D subscale sampler. Theneach sampler filters out all but its subscale points. FIG. 12illustrates two possible multi-sampler designs.

1.2 Digital Filter Bank

The heart of the D Digital Filter design is a bank of ordinary,commercially available IIR or FIR digital filters, as shown in FIG. 2.Each filter in the bank is associated with a subscale of the underlyingD Scale used for sampling. That is, each filter independently samplesthe incoming data at the rate of its associated D subscale.

A 1 dimensional (1D) D digital filter is characterized by a vector ofvectors of complex valued parameters or tap weights. This vector ofvectors corresponds to the individual D subscales in the underlying DScale used for the non-uniform sampling. It is realized by a bank of theusual digital filters, with each characterized by one vector of tapweights from the vector of vectors. The weights are the coefficients ofthe discrete impulse response.

1.3 Sample Router

Additionally, a router component directs each new sample to a filter inthe filter bank. This generalizes the usual 1D digital filter that ischaracterized by a vector of weights. FIG. 13 illustrates thearchitecture. FIG. 14 provides a close up of the filter bank router.This component first maps each incoming point into a D Scale point, ifthe point is not already a D Scale point. It then routes the mappedpoint to that filter in the D filter bank corresponding to the Dsubscale containing the mapped point.

1.4 Multi-Stream Consolidator

This component consolidates the multiple D subscale uniform samplestreams into one stream consisting of interleaved samples from themultiple D subscale streams. The net result is one stream of samplesflowing at a non-uniform rate.

1.5 Phase Shifter

This component changes the phase of processed D subscale point streams(or signals). It performs this task by delaying the passage of each Dsubscale point by a predetermined value. In effect, the phase shiftertransforms a point on one subscale into a point on another subscale oreven onto the same subscale, but at a different frequency. A phaseshifter might be useful for applications involving positive and negativefeedback loops.

1.6 Sampling Feedback Loop Controller Component

The feedback loop can selectively feedback control messages or signalsto the D Scale sampler at the input stream. These control signalseffectively modulate the input stream that actually passes through thechannel being filtered. This mechanism adds a new input samplingstrategy, in addition to random and trend sampling. FIGS. 6 and 9illustrated where this component can be integrated into the D filterdesign.

1.7 D Arithmetic Tables

The D Arithmetic is detailed in patent application Ser. No. 09/568,368.Actual floating point measurements or sampled values are mapped to DScale points. Arithmetic operations are performed on these points bymapping the floating point result of the arithmetic operation on eachpair of operands into another D Scale point. Scaling is properlypreserved as well when the D Arithmetic tables are created. A table canbe a physical region in memory of a computer host. Or it can be alogical table of references or pointers to a network of physical tablesdistributed across multiple computers.

The D Arithmetic tables can be readily integrated into a filteringsystem to perform the calculations, which become iterative mappingsacross D Arithmetic table mappings. The D Arithmetic is portable and itdoes not have floating point operation overflow, underflow ordivide-by-zero conditions. Further, it has a self-stabilizing propertythat effectively bounds round-off errors to within a tolerance value.This tolerance depends on the minimum possible distance between alladjacent points in the D Scale that underlies the D Arithmetic tables.

Therefore, use of the D Arithmetic in filter design can diminish, if notmitigate entirely, the corrosive effects of floating point calculationsin filter design.

2. Rationale for D Digital Filter Bank Design

2.1 Decomposition of Convolution Integral

The motivation, indeed, the underlying mathematical justification, forthis design, is presented followed by details on determining thefilter's tap weights.

As mentioned previously in the background section and as shown in FIGS.2 and 3, the vector of weights in an industry standard digital filterderives from the standard discrete approximation of the convolutionintegral of a signal and the impulse response of a channel.

These N complex coefficients are the tap weights in a transversaldigital filter. The rationale for this linkage of coefficients to filtertaps is that the vector of tap weights approximates the impulse responseof a linear system.

However, the D Scale can be applied to calculate the exact convolutionof a complex function. Consider the same linear system illustrated inFIG. 1. Let the D Scale be comprised of N subscales {p₁, . . . p_(N)}.

Let δ_(i,j) denote the distance of point p_(j) in D subscale i to thenext D Scale point. Recall that the sampling is uniform at a rate of1/p_(i) for each subscale i. Therefore, for points on subscale i:

t=n/p _(i) where n=1, . . . p_(i)−1

(Note that for the smallest scale, p_(min) n starts at 0. That is, theconvention is that the smallest scale contains the 0 D Scale point).

Further, the convolution integral can be decomposed using Riemannsummation over the D subscales, identically to the procedure used toderive an exact decomposition of the Fourier Transform in patentapplication Ser. No. 09/326,084. The continuous convolution integral iscalculated by determining its value separately at points on each Dsubscales.

On each subscale p_(j), the integration variable τ extends over uniformintervals 1/p_(j). Therefore on each subscale:

τ=k/p _(j)

Also, on each subscale, the non-uniform interval associated with eachuniformly distributed D Scale point leads to:

$\begin{matrix}\left. {\tau}\rightarrow\delta_{j,k} \right. \\{{y(t)} = {{\int_{- \infty}^{+ \infty}{{h\left( {t - \tau} \right)}{u(\tau)}{\tau}}} = {{\int{{h(\tau)}{u\left( {t - \tau} \right)}{\tau}}} \approx}}}\end{matrix}$

Decomposition of Convolution into D Subscale Filters

$\begin{matrix}{{y\left( {n/p_{i}} \right)} = {\sum\limits_{j \in {\{{p_{1},{\ldots \; p_{N}}}\}}}{\sum\limits_{k = 0}^{k = {p_{j} - 1}}{{h\left( {k/p_{j}} \right)}{u\left( {{n/p_{i}} - {k/p_{j}}} \right)}\delta_{j,k}}}}} \\{{n = 1},{{\ldots \mspace{14mu} p_{i}} - 1}}\end{matrix}$

This expression is then used to calculate the continuous convolutionintegral across the range of interest, by assembling the calculatedconvolution values along points on each D subscale. Again, thisprocedure is identical to that described in patent application patentapplication Ser. No. 09/326,084 to calculate an exact decomposition ofthe Fourier Transform of a transient signal. The approximation becomesmore exact as more D subscales are used, which increases the resolutionof the D Scale.

2.2 Internal Architecture of D Filter Bank

Note that this discrete approximation to the continuous convolutionintegral contains N of the usual discrete impulse response convolutionseries, one for each D subscale. This is the key observation toestablishing the link between the coefficients of these series and the Dfilter bank tap weights. For consider the following slight rearrangementof the terms in this equation: Let a weight be denoted by:

w* _(j,k) =h(k/p _(j))δ_(j,k)

where the asterisk denotes complex conjugation, which assumes that thetap inputs, and therefore the tap weights are all complex valued. Then:

${y\left( {n/p_{i}} \right)} = {\sum\limits_{j \in {\{{p_{1},{\ldots \; p_{N}}}\}}}{\sum\limits_{k = 0}^{k = {p_{j} - 1}}{w_{j,k}^{\star}{u\left( {{n/p_{i}} - {k/p_{j}}} \right)}}}}$

The equation can now be interpreted in the language of filter weights.For each internal summation convolves the finite duration impulseresponse of one filter, represented by the vector w*_(j,k) for Dsubscale j, with the input u(n/p_(i)) from D subscale i, to produce anoutput contribution y(n/p_(i)) for a point on D subscale i. FIG. 15illustrates the internal structure of an FIR D filter bank that realizesthis equation. Note that the conventional unit delay box z⁻¹ has anadditional subscript associated with the D subscale. z₀ ⁻¹. The intentis to highlight the variable “unit” interval as it depends on the Dsubscale.

The D filter bank of IIR filters follows the same pattern whereby the DScale δ_(j,k) values are multiplied with the original weights and a setof N conventional IIR filters are tied to a D Scale sampler and router.

2.3 D Signal

2.3.1 Motivation

The preceding equations assumed that the impulse response of a systemwas already known. There is a then a need to actually determine ormeasure the impulse response of a system (or channel). The standardtheoretical approach is to discretize the Dirac functional as follows:

${h(t)} = {{\int_{- \infty}^{+ \infty}{{h(\tau)}{\delta \left( {t - \tau} \right)}{\tau}}} \approx {\sum\limits_{n = {- \infty}}^{+ \infty}{{h\left( {n\; {\Delta\tau}} \right)}{\delta \left( {t - {n\; {\Delta\tau}}} \right)}{\Delta\tau}}}}$

The effect of any input x(t) to the linear system is then determinedusing the system's linearity property as:

${y(t)} = {{\int_{- \infty}^{+ \infty}{{h(\tau)}{x\left( {t - \tau} \right)}{\tau}}} \approx {\sum\limits_{n = {- \infty}}^{+ \infty}{{h\left( {n\; {\Delta\tau}} \right)}{x\left( {t - {n\; {\Delta\tau}}} \right)}{\Delta\tau}}}}$

The approximation of the continuous Dirac impulse response is a spike ofunit amplitude. Note that phase is ignored in this procedure.

The motivation for the D Signal is the set of equations in patentapplication Ser. No. 09/326,084 on an exact decomposition of themagnitude and phase of the Fourier transform of a signal usingnon-uniform D Scale sampling. The equation for the magnitudedecomposition is:

${{F\left( {2\pi \; {m/{Tp}_{i}}} \right)}} = {\sum\limits_{j \in {\{{p_{1},{\ldots \mspace{14mu} p_{N}}}\}}}{\sum\limits_{k = 0}^{k = {{{Tp}_{k} \star p_{k}} - 1}}{\delta_{{p\; k},{{n/p}\; k}}{h\left( {n/p_{i}} \right)}^{{- 2}{\pi {({{Tp}_{k}/T_{pi}})}}n\; {m/{({T_{p\; k} \star p_{k}})}}}}}}$

The decomposition consists of a double sum of terms with each term ofthe form:

δ_(pk,n/pk)h(n/p_(i))e^(−2π(Tp) ^(k) ^(/T) ^(pi) ^()nm/(T) ^(pk) *^(p)^(k) ⁾

The inverse Fourier Transform decomposition is used to determine theshape of a signal derived from a specified frequency spectrum.

2.3.2 Relation to Impulse Function

These decomposition equations can be interpreted as a formula formeasuring the impact of a channel or system on linear combinations ofcertain atomic or elemental signals. An elemental signal is one term inthe decomposition, whether forward or inverse.

Let a D Signal be one of the terms in the decomposition above. Then any(reasonably well behaved) signal can be decomposed into a set of Dsignals. The D signal, then realizes a D Scale point. It ischaracterized by a predominant monotonic frequency at the rate of thesubscale's generator. If prime numbers are used, then the dominantfrequency is 1/p_(k) where k is the k^(th) subscale in the D Scale used.In practice, there will be a spread, quantified by a variance, in thefrequency spectrum of the D Signal. This variance must be within themaximum distance between adjacent points of the D Scale, as shown inFIG. 16. This insures that each D Scale point's corresponding D Signalis distinct.

Thus, if a set of signals is desired with a flat frequency spectrum in aspecified range, then the atomic terms can set their middle factor h( )to 1. Further, if the phase is flat, then the exponential termsimplifies considerably. The net effect of these constraints is a signalwith a flat frequency and phase response. Indeed, by definition, thatsignal would realize the impulse in the defined frequency band. The DSignal approximation of an impulse response can be tailored to thefrequency domain of interest, as in the microwave or optical frequencyranges.

2.3.3 Possible Uses of D Signals

One possible use of the D Signal is as CDMA signal components. Anotheruse is to empirically measure the impulse response of a possiblynon-linear system, discussed in more depth later.

2.4 Nonlinear Systems Addressed by D Filter Design Methodology

Nonlinear systems exhibit a wide range of behavior. Consider that oneend of the spectrum of behavior is almost linear (i.e, superposition andtime invariance apply) while at the other end is extreme behavior suchas expected at the rim of a black hole. The scope of nonlinear systemsconsidered by the D filter design is bounded by the part of the spectrumto include non-time varying systems that exhibit deterministic,multi-frequency generation per monotonic or monochromatic frequencyinput.

The following methodology can be used to quantify such nonlinearities.First, a D Scale is created for the application domain considered, suchas microwave, optical or millimeter regions. Then each unit amplitude DSignal associated with the D Scale is applied at the input. Because ofthe system's nonlinearities, outputs at multiple frequencies can beexpected, as shown in FIG. 17. These output signals will have variedamplitudes and phases. And their frequencies can be mapped to D Scalefrequency points. Again, the nonlinearities are assumed to be timeinvariant to enable reproducible behavior. The net result of theseefforts is that the output spectrum can be expressed as follows,illustrated in FIG. 18A:

Y _(i,j)(ω)=H _(i,j)(ω)*X _(j)(ω)

Where:

Y_(ij) is the component of output D Signal associated with subscale idue to an input D Signal from D subscale j.And the following results after collecting contributions from all inputD Signal frequencies, as shown in FIG. 18B:

$\underset{j = 0}{\overset{N - 1}{Y_{i}(\omega)}} = {{\sum\underset{j = 0}{\overset{N - 1}{Y_{i,j}(\omega)}}} = {\sum{{H_{i,j}(\omega)} \star {X_{j}(\omega)}}}}$

N=number of D subscales0≦i<N−1

Where:

-   -   H_(i,j)→ Transfer function for input D signal associated with        subscale j and output D Signal associated with subscale i.    -   X_(j)→ Input D signal associated with subscale j.    -   Y_(i)→ Output D Signal associated with subscale i.        Or in matrix notation,

Y=H*X

Where:

H is an N×N matrix each of whose elements is a transfer functionH_(i,j)(ω)

X is a 1×N vector of input D Signals

Alternately, this can be expressed in the time domain. In that case,there are many impulse response functions to consider, one for each DSignal.

Note again, that for this general case contains the usual relationshipsbetween input and outputs, when the system is linear. For then, there isonly one output frequency per input D Signal and the transfer functionmatrix only contains diagonal elements That is,

Y _(i)(ω)=H _(i,i)(ω)*X _(i)(ω)

Or as usually expressed in linear system theory:

Y(ω)=H(ω)*X(ω)

Nonlinearities might also occur such that the change in output amplitudeis not proportional to input amplitude changes. In that case, the sametest methodology can be applied, but varying input amplitudes on eachtrial. That is, after a set of output signals is collected for a rangeof D Signal inputs at a set amplitude, a new set of D Signal inputs isapplied, but with the amplitudes incremented. This procedure isillustrated in FIG. 19.

2.5 Method for Empirically Measuring Impulse Response of a Linear orNonlinear System

The previous discussion can be summarized through the followingprocedural description. A flow chart illustrating the methodology ispresented in FIG. 20.

2.5.1 Construct a D Scale for the Problem at Hand:

Decide on which subscales to include in a D Scale based on a-priorknowledge of the application domain.

2.5.2 Generate D Signals

Generate the D Signals corresponding to the D Scale points constructedin the previous step. Just use the steady state part of the signal.

2.5.3 Apply D Signals

Apply the D Signals to the channel of interest, to determine the impulseresponse.

2.5.4 Measure the Responses for Each Signal Application.

There is a matrix of input signals. Note that in the D Scale range,nonlinearities can create multiple output frequency results per inputfrequency. Therefore each output is the composition of multiplecontributions from D Signals.

2.5.5 Calculate the Transfer Function Matrix as Follows:

H _(i,j)(ω)=Y _(i,j)(ω)/X _(j)(ω)

2.5.6 Calculate the Impulse Responses by Taking the Fourier Transform ofthe Calculated Transfer Functions.

The D Scale decomposition should be applied, rather than the standard,coarse, DFT approximation. This effort will yield a vector of vectors ofpossibly complex coefficients.

3. Applications of D Filter Design

3.1 Linear Static Filter Design Algorithms

There are many well known and applied categories of linear filter designalgorithms, and many implementation approaches within each category. Forexample, one approach is frequency domain sampling, while anotherapproach, typically used for IIR designs, relies on complex planemappings.

State-of-the-art digital filter design approaches typically perform twogeneral steps. First is a determination of the digital filterarchitecture to use, based on the application requirements and resourceconstraints. The two basic architectures are Finite-Impulse-Response(FIR) or Infinite-Impulse-Response (IIR) filters. Second, a vector ofpossibly complex weights or transverse filter taps is calculated.

Each existing technique is readily generalized to the D filter design.For the D filter design is a direct generalization of the standardvector of tap weights. For example, the frequency sampling approach isgeneralized to use non-uniform D Scale sampling of a frequency spectrum,rather than the uniformly spaced samples usually applied. The same stepsare then followed for the set of uniformly sampled D subscales. Thissame reasoning can be applied to each current digital filtering designapproach.

Standard filter design entails the determination of the vector of tapweights of a filter. The three most commonly used FIR design approachesare, along with their D filter generalizations:

1. Frequency Sampling

The desired frequency response is sampled uniformly at a sufficientlyhigh rate. Then the discrete inverse DFT of the samples is calculated.

This technique is readily generalized for the frequency sampling can beperformed non-uniformly using the D Scale. Then the inverse Fouriertransform can be more accurately calculated using the D Scale samples.

2. Windowing

This technique entails multiplying the desired frequency response by a“windowing” function. The purpose is to reduce the impact of the Gibbsphenomenon when approximating the impulse response using the discreteDFT.

This technique is also readily generalized. For the inverse calculationcan be performed more accurately allowing for more varied windowingfunctions.

3. Min-Max Optimization or Least-Mean-Square Error

This technique minimizes a mean-squared error criteria by solving aderived set of linear equations. Many algorithms are in use includingLevison-Durbin, Parks-McClellan and Remez exchange, to note some of themost widely used ones.

This technique is also readily generalized for D filter design for theerror measure would contain extra sets of terms, each set correspondingto a D subscale. This will be discussed further in the followingsections, as well as in the section that discusses the derivation of thegeneralized Weiner-Hopf equation. Optimizing that error measure thenentails the concurrent solution of multiple sets of linear equations,with each set associated with a D subscale. An attraction is that eachset of linear equations can be of an order that is possiblysignificantly lower than that of the standard technique. This is due tothe reduced number of sample points (or sampling rate) for each Dsubscale compared to the far higher rate needed by the usual technique.

The most commonly used IIR filter design technique is:

1. Bilinear Transformation.

In this general approach, a transfer function is represented as afraction containing poles and zeros, for the most general IIR case. Thegeneralized D filter transfer function can also be expressed as either aseries of such fractions, or combined into one fraction. Then the samestandard techniques could be applied.

Several other techniques, less commonly known or used, are available todesign filters using both magnitude and phase specifications.

3.2 Nonlinear Filter Design Algorithms

The same enhancements to the standard design techniques apply to thenonlinear systems described in Section 2.4. For this class of nonlinearsystems can be characterized by an extension of the sum of sum of tapweight vectors described previously. Therefore the extra vectors in thenonlinear systems are accommodated by simply extending the descriptionsin the previous Section 3.1.

3.3 Filters Derived from the Weiner-Hopf Equation

3.3.1 Conventional Weiner-Hopf Equation

The Weiner-Hopf equation underlies the design of adaptive filters aswell as of recursive, predictive, Kalman filter types. It is derivedfrom a linear model that relies on uniformly sampled input and outputvalues, shown in FIG. 21. An equation is derived by seeking the extremepoint of a least-mean-square expression to yield, using the nomenclaturefrom “Adaptive Filter Theory” by S. Haykin, 3^(rd) edition,

Rω=P

where:

R=E[u(n)u*(n)]]→M×M autocorrelation matrix of input

u(n)→M×1 tap input vector [u(n), . . . u(u−M+1)]^(T)

ω=[w₀, . . . w_(M−1)]^(T)→M×1 vector of transversal filter's tap weightthat optimize an error criteria.

P=E[u(n)d*(n)]→M×1 cross correlation matrix of input-output

The solution is a vector of possibly complex valued tap weights in thetransverse filter. This equation is solved iteratively by descending anM+1 dimensional error performance surface that is shaped as a bowl. Ithas M degrees of freedom represented by the filter's tap weights. Theoptimal solution is at the bottom of the bowl, as shown in FIG. 22. Ateach step:

(Updated value of)=(Last value of)+(Learning rate)*(Tap inputvector)(tap weight)(tap weight)(parameter)(error signal)

This recursively determined solution is realized by dynamically changingtap weights in an adaptive filter.

3.3.2 Generalized Weiner-Hopf Equation for Linear Systems

A generalized Weiner-Hopf equation for linear systems is derived from amathematical model that uses non-uniformly sampled input and outputsamples, as well as the D Signal based impulse response, described inSection XXX. The standard Weiner-Hopf equation and model is shown as aspecial case when uniform sampling is applied. The generalized equationforms the foundation for D filter design.

The generalized Weiner-Hopf equation is derived using the samemethodology as that for the standard equation. FIG. 23 presents themodel used by the derivation. The derivation relies on the observationthat the optimal value of a set of series of terms is located byoptimizing each series of terms in the equation. In this case, the setof series occurs because each is associated with one of the D subscales.What results is a set of smaller dimensional conventional Weiner-Hopfequations. The set of such equations can be organized into a supermatrix Weiner-Hopf equation to yield:

Generalized Weiner-Hopf Equation

R ω= P

where:

-   -   R=E[u_(pi)(n)u_(pj)*(n)]]→Matrix of N×N autocorrelation matrices        of D subscale inputs. Each element is an autocorrelation matrix        whose elements are correlations between input points on D        subscale i and j, for all D subscales. N is the number of        subscales in the D Scale.    -   u_(pi)(n)→M_(i)×1 tap input vector [u_(pi)(n) . . .        u_(pi)(u−M_(i)+1)]^(T) where each sample is in D subscale i.        M_(i) is the number of samples that depends on the D subscale i.    -   ω=[W₀, . . . W_(N−1)]^(T)→Vector of N×1 vectors W_(i)=[w_(i0), .        . . w_(iMi−1)]^(T) Each vector contains the tap weights for the        transversal filter in a D filter bank associated with D subscale        i.    -   P=E[u_(i)(n)d_(j)*(n)]]→Matrix of N×1 cross correlation matrices        of D subscale input-outputs. Each element is a cross correlation        matrix whose elements are correlations between input and output        points from D subscales i and j respectively.    -   i and j→0<=i,j<N where N is the number of D subscales.

The solution is a vector of vectors of possibly complex valued tapweights in each D filter bank transverse filter. This equation is solvedby simultaneously solving each of its standard Weiner-Hopf matrixequation components. Each component is associated with a pair of Dsubscales, i and j. As with the standard equation, each componentequation is solved iteratively by descending an M_(i)+1 dimensionalerror performance surface that is shaped as a bowl, where M_(i) is thenumber of samples in D subscale i. It has M_(i) degrees of freedomrepresented by the filter's tap weights. The optimal solution for eachequation is at the bottom of each bowl, as shown in FIG. 24. Again, aswith the standard case, at each step, for each set of tap weights in theD filter bank:

(Updated value of)=(Last value of)+(Learning rate)*(Tap inputvector)(tap weight)(tap weight)(parameter)(error signal)

This recursively determined solution is realized by dynamically changingtap weights in each adaptive filter in the D filter bank.

3.3.3 Discussion of Generalized Weiner-Hopf Equation

The standard equation results when the special case of uniform samplingis applied. For then, the D Scale has only one subscale and the matrixof vector of vectors collapses to the standard form. Indeed, theunderlying model then reverts back from that shown in FIG. 23 to thestandard model shown in FIG. 21.

Another point worth noting is that the generalized equation has astructure that is reminiscent of the fractional adaptive equalizer. Buteven as the syntax describing both systems is similar, the semantics ofthe two differ considerably. One difference is that the fractionalequalizer uses a higher data rate whereas the D Filter is designedaround lower rate D Scale signals.

The interpretation of the solution of this equation is a generalizationof the interpretation of the usual equation. The search for a solutionof this equation can be interpreted as the simultaneous descent downmultiple bowls of possibly significant lower dimensions. This isillustrated in FIG. 21. Each extrema of each bowl corresponds to avector of possibly complex tap weights of a filter in the D filter bankassociated with a D subscale. The full solution of the generalizedWeiner-Hopf equation then consists of the vector of vectors, one foreach of the D subscales.

3.3.4 Generalized Weiner-Hopf Equation for Nonlinear Systems

Nonlinearities of the type described in Section XXX can be readilyaccommodated by the generalized Weiner-Hopf equation. For thenonlinearities in a channel add additional summations of terms to theimpulse response characterization. And these added summations manifestthemselves by additional elemental, conventional, Weiner-Hopf equationsthat must be solved simultaneously. The reasoning is identical to thatwhich yielded the original generalized Weiner-Hopf equation for linearsystems, as discussed previously.

The interpretation of the solution for a nonlinear system is of anadditional number of bowls, but each a higher dimension, determined bythe extent of the nonlinearities.

3.4 Benefits of Generalized Weiner-Hopf Based Adaptive Filters

The generalized Weiner-Hopf equation underlies the enhanced design ofadaptive filters. Its recursive solution is realized by a set ofadaptive filters whose tap weights are varied according to the D Scalenon-uniformly sampled inputs. Adaptive filters that are based on Dfilter design have the following advantages over comparableimplementations that rely the standard model.

-   a. Reduced number of taps because each D subscale has far fewer    weights than the uniformly sampled case.-   b. More sensitive to variations in input signal and able to track    the variations and volatility. This falls out of the more accurate    characterization of a channel.-   c. Non-calculating or pre-wired D Arithmetic to avoid the    deleterious effects of round-off errors in real-time filter    calculations.-   d. Possibly reduced number of iterations for convergence or possibly    the time to convergence. This is because the descent down a deep    high of order N dimensional bowl is replaced by simultaneous    descents down several lower dimensional order and shallow bowls.-   e. Readily accommodates certain nonlinearities without linearizing    away the problem or approximating it.-   f. It enables IIR solutions due to enhanced stability.-   g. It enables the use of multiple sets of transversal filters, each    with potentially far fewer number of taps than the conventional    filter. Fewer taps along the end-to-end path of a sample's traversal    is generally advantageous.

3.5 Computing Auto Correlation and Cross-Correlation

As with the industry standard digital filter, practical implementationof digital filters requires the accurate computation of the inputsignal's autocorrelation as well as of the cross correlation of theinput and output signals. The same D Scale decomposition technique usedto derive an exact numeric expression for the Fourier (or other kernel)transform can also provide accurate numeric calculations ofautocorrelation and cross-correlation functions of possibly complexvalued signals. Further, the derivation is almost identical to thatpresented previously to compute the convolution.

Consider the continuous autocorrelation integral for an input signal:

R_(uu)(t) = ∫_(−∞)^(+∞)u(t)u(t + τ)τ

Let the underlying D Scale used to sample a signal be comprised of Nsubscales {p₁, . . . p_(N)}.

Let δ_(i,j) denote the distance of point p_(j) in D subscale i to thenext D Scale point.

Recall that the sampling is uniform at a rate of 1/p_(i) for eachsubscale i. Therefore, for points on subscale i,

t=n/p _(i) where n=1, . . . p_(i)−1

Further, the correlation integral can be decomposed using Riemannsummation over the D subscales. The continuous correlation integral iscomputed by determining its value separately at points on each Dsubscales.

On each subscale p_(j,) the integration variable τ extends over uniformintervals 1/p_(j). Therefore on each subscale:

τ=k/p _(j)

Also, on each subscale, the non-uniform interval associated with eachuniformly distributed D Scale point leads to:

dτ→δ_(j,k)

Finally, after substituting these mappings into the Riemann summationthat computes the original integral, the following equation emerges:

Computation of Autocorrelation using the D Scale:

${R_{uu}\left( {n/p_{i}} \right)} = {\sum\limits_{j \in {\{{p_{1},{\ldots \mspace{14mu} p_{N}}}\}}}{\sum\limits_{k = 0}^{k = {p_{j} - 1}}{{u\left( {n/p_{i}} \right)}{u\left( {{n/p_{i}} + {k/p_{j}}} \right)}\delta_{j,k}}}}$0 ⇐ i < N

Note that, as in the computations of the Fourier transform andconvolution, the final result is generated as by concatenating thepartial results for each D subscale.

Similarly, consider the continuous cross-correlation integral forinput-output signals:

R_(uy)(t) = ∫_(−∞)^(+∞)u(t)y(t + τ)τ

Computation of Cross-correlation using the D Scale:

${R_{uy}\left( {n/p_{i}} \right)} = {\sum\limits_{j \in {\{{p_{1},{\ldots \mspace{14mu} p_{N}}}\}}}{\sum\limits_{k = 0}^{k = {p_{j} - 1}}{{u\left( {n/p_{i}} \right)}{y\left( {{n/p_{i}} + {k/p_{j}}} \right)}\delta_{j,k}}}}$0 ⇐ i < N

Note again, that, as in the computations of the Fourier transform andconvolution, the final result is generated as by concatenating thepartial results for each D subscale.

The autocorrelation functions can also be computed via the indirectapproach of computing the spectrum of the signal and relying on theWeiner-Khintchine Theorem. The D Scale based decomposition of thespectrum would then be applied.

3.6 Generalized Weiner (Matched) Filter

3.6.1 Conventional Weiner Filter

The Weiner filter is also based on the Weiner-Hopf equation although itis typically presented as a matched filter. It relies on a-prioriknowledge of to correlate the combined input signal and noise againstthe expected signal. This why it is often referred to as a matchedfilter, for it aims to match the input to a specific signal. FIG. 23shows the standard continuous Weiner filter. FIG. 24 shows the usualdiscrete realization of the continuous formulation. There, onecorrelator is applied to uniformly sampled inputs to compute:

${r_{uy}(I)} = {{1/N}{\sum\limits_{n = 0}^{N - 1}{{u(n)}{y\left( {n + I} \right)}}}}$

where 1/N is the sample rate.

3.6.2 Generalized Weiner Filter

The generalized Weiner filter design uses a D Scale data samplerfollowed by a router to a bank of industry standard Weiner matchedfilters. The signal to be matched at each Weiner filter in the bank issampled at the rate of its corresponding D subscale. This is illustratedin FIG. 25.

Note that the design is a subset of the most general D filter designarchitecture shown in FIG. 6. Note also from that figure, that there areN vectors of weights. Each weight corresponds to a coefficient in theexact autocorrelation expansion of the continuous autocorrelationfunction of the input signal, however transient or non-stationary.

3.7 Generalized Kalman Filter

The Kalman filter is also based on the Weiner-Hopf equation. But it ismost often formulated in the state variable formalism that make theunderlying foundation less obvious. The standard Kalman filter isderived from a pair of process (or state) and measurement equations. Asignal flow graph representation of a linear discrete-time dynamicsystem is shown in FIG. 26.

That figure illustrates the following equations:Process (or state) Equation:

x(n+1)=F(n+1,n)*x(n)+v ₁(n)

-   x(n)→ Unknown state vector of dimension M at time n unit intervals.-   F(n+1,n)→ Known M×M state transition matrix relating the state of    the system at times n+1 and n unit intervals.-   v₁(n)→ M×1 vector representing process noise at time n unit    intervals

Measurement Equation:

y(n)=C(n)*x(n)+v ₂(n)

-   C(n)→ Known N×M measurement matrix.-   v₂(n)→ N×1 vector representing measurement noise at n unit    intervals.-   e⁻ _(k)=x(n)−x⁻(n)→ Estimation error between state vector and prior    estimate or best estimate about the process prior to time n.-   P_(k)=E[e⁻ _(k)e⁻ _(k) ^(T)]→ Predicted state error covariance or    correlation matrix.

The Kalman filtering problem is to use the entire observed data set,consisting of vectors y(1), y(2), . . . y(n), to find for each n≧1, theminimum mean-square estimates of the components of the state x(i).

The generalized Kalman filter is derived using the same analytical stepsto derive the standard Kalman filter. The derivation is based on thegeneralized signal flow graph representation shown in FIG. 27. Thatfigure illustrates the flow realization of the following equations:

Process (or state) Equation:

x _(j)(l)=F _(i,j)(l,k)*x _(i)(k)+v _(1i)(k)

-   x_(i)(k)→ Unknown state vector of dimension M from D subscale i at    time k intervals on that D subscale.-   x_(j)(l)→ Unknown state vector of dimension M from D subscale j at    time l intervals on that D subscale. This is the very next point on    the D Scale following x_(i)(k).-   F_(i,j)(l,k)→ Known M×M state transition matrix relating the state    of the system at adjacent points on the D Scale, that at I intervals    on D subscale j and the previous point at interval k on D subscale    k.-   v_(1i)(k)→ M×1 vector representing process noise at n unit    intervals.

Measurement Equation:

y _(i)(n)=C _(i)(n)*x _(i)(n)+v _(2i)(n)

-   C_(i)(k)→ Known N×M measurement matrix at a point on the D Scale    that is at the k-th interval on D subscale i.-   v_(2i)(k)→ N×1 vector representing measurement noise at a point on    the D Scale that is at the k-th interval on D subscale i.-   e⁻ _(k,i)=x_(i)(k)−x⁻ _(i)(k)→ Estimation error between state vector    and prior estimate or best estimate about the process prior to this    D Scale point.-   P_(k,) _(i) =E[e⁻ _(k,i)e⁻ _(k,i) ^(T)]→ Predicted state error    covariance or correlation matrix.

The generalized Kalman filtering problem is to use the entire observeddata set, consisting of vectors at times on the D Scale y_(i)(1), . . .y_(j)(l), to find, for each time along the D Scale (i.e., each partialtime series), the minimum mean-square estimates of the components of thestate x_(i)(k) on each D subscale.

The only difference in the analysis, is that each new sample is mappedto a D Scale point that is associated with a unique D subscale. Thecorrelation matrix element in the super-matrix of the generalizedWeiner-Hopf equation is identified from its D subscale. The correlationmatrix element corresponding to that D subscale must be used for thatcycle in the Kalman filter's recursive algorithm. Therefore the onlyaddition to the Kalman iterative procedure is to identify whichcorrelation matrix in the super-matrix to use on each iteration.

The resulting modification to the standard Kalman Filter is shown in theflow chart in FIG. 28, with the new steps highlighted. Again, the extrastep is to fetch the corresponding correlation matrix element for eachnew sample point.

Note again that the usual Kalman filter algorithm appears for thespecial case of uniform sampling. Also note that nonlinearities aretransparently factored into the calculations because their effects arealready embedded into the off-diagonal weights and their correspondingcorrelation matrix elements.

1-41. (canceled)
 42. A method of modulating a nonuniformly sampled datastream to the digital filter using a controller, comprising usingcontrol signals based on the output, applying the control signals toinput samplers, driven by its local volatility calculation.
 43. Themethod of modulating a nonuniformly sampled data stream to the digitalfilter using a controller of claim 42 further comprising calculating thelocal volatility of adjacent nonuniform samples from the D Scalesampler, via the slope, convexity and/or n-th local derivative using Nof the last nonuniform samples, where N depends on the sensitivity ororder of the derivative desired.
 44. The method of modulating anonuniformly sampled data stream to the digital filter using acontroller of claim 42 wherein the controller applies control signals tosuppress the D Scale's subscales' uniform samplers for specific samplingresolutions or frequencies.
 45. The method of modulating a nonuniformlysampled data stream to the digital filter using a controller of claim 42wherein the controller applies control signals to trigger the S Scale'ssubscales' uniform samplers for specific sampling resolutions orfrequencies.